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Abstract 



o 
O 

^— «j , Phase reduction is a commonly used techinque for analyzing stable oscillators, par- 

' ticularly in studies concerning synchronization and phase lock of a network of oscilla- 

tors. In a widely used numerical approach for obtaining phase reduction of a single 
oscillator, one needs to obtain the gradient of the phase function, which essentially 
' provides a linear approximation of isochrons. In this paper, we extend the method for 

, obtaining partial derivatives of the phase function to arbitrary order, providing higher 

• ■ order approximations of isochrons. In particular, our method in order 2 can be applied 

^ ' to the study of dynamics of a stable oscillator subjected to stochastic perturbations, a 

topic that will be discussed in a future paper. We use the Stuart-Landau oscillator to 
illustrate the method in order 2. 

> ■ 1 Introduction and statement of main results 

O 

On . Weak perturbations of limit cycle oscillators are of great interest in a variety of fields in 

physics, chemistry, engineering, and quantitative biology, whenever the system under study 
displays stable oscillations. A powerful theoretical approach in the analysis of weakly per- 
turbed limit cycles, particularly in relation to synchronization and phase lock of a network 
of oscillators, is to reduce the description of the system to a single "phase" variable. This 
phase reduction procedure is the focus of the present paper. 

In this introduction, we begin by describing a standard numerical approach for obtaining 
first order phase reduction, due to I.G. Malkin [TlJ|T2]. We then explain our higher order 
, method, which is summarized in Theorem 11.11 and illustrate it in order 2 using the well- 

^ ' known Stuart-Landau oscillator as an example. Proofs are are given in the subsequent 

sections. Our focus is on the theoretical underpinnings of the method rather than on the 
details of numerical implementation, but we provide a numerical example to illustrate the 
approach. 

The following set-up is assumed to hold throughout. Let _F be a smooth (i.e., continu- 
ously differentiable to all orders) vector field on n-dimensional Euclidean space M" . The flow 
line of an initial point x will be denoted (j)t{x), or simply x{t), for t G R. Let 6 be a stable 
(hyperbolic) limit cycle of the differential equation x = F{x) having period T > 0. We 
write to = 1/T for the reciprocal of the period. The stability condition means that for any 
given point x on 6, there exists a linear n — 1-dimensional subspace W{x) in M" transverse 
to F{x) such that vectors in W{x) contract exponentially under positive iterations of the 
differential {d<t)T)x of the flow map at x. That is, |((i0T)"u| < CA" |v| for positive constants 
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C, A < 1, positive integer n, and all v in W{x). As we are only concerned with the system 
near 6, there is no loss of generality in assuming that is a complete vector field, so that 
the flow lines (j>t{x) are defined for alH G M and all x € M". 

1.1 The phase function 

We briefly review a few well-known facts about the dynamics near a stable limit cycle for 
the purpose of setting up notation. 

Due to normal hyperbolicity on C, it is known [6] that a neighborhood of the limit cycle 
is continuously foliated by contracting manifolds, 'W(x), for each x G C, where the W(x) are 
smooth submanifolds tangent to W{x) at x and diffeomorphic to an open disc. The foliation 
is invariant, i.e., (pt maps W{x) into W{(j)t{x)) for all t > and all x in C. Furthermore, the 
foliation is smooth since C is an orbit of a smooth flow. In fact, we can use maps of the 
form H{z,t) = (pt{z), where z £ 'W{x) and t G (—a, a), to produce smooth foliation charts 
around x. 

Let be a function defined on a sufficiently small neighborhood of C, whose level sets 
are the local contracting manifolds and such that Q{x{t)) = tot modulo integer translations. 
I.e., 8 takes values in R/Z. We often regard Q simply as a function into R, keeping in mind 
that, in this case, 8(x(i + nT)) = Q{x{t)) + n. Derivatives of O, on the other hand, are 
single-valued functions into M. 

We refer to as the phase Junction of the oscillator. It is, thus, a smooth function on 
some neighborhood of the limit cycle such that 0(j/) = 8(x) iff \4)t{y) — (j)t{x)\ i) ast ^ oo 
whenever a; S C and y lies in a sufficiently small neighborhood of C. 

The level sets of 8 are also called isochrons [18j. The existence of isochrons was first 
proved in [4] . Since the component functions of the gradient of 8 give the change in phase 
due to a small change in the respective position coordinates [1] , the graphs of those functions 
are often referred to as phase response curves. This notion is widely used in theoretical and 
experimental neuroscience. For example, some theoretical studies have examined how the 
shape of the phase response curves affect the sychronization dynamics of coupled oscilla- 
tors [3 [2] . Other studies have investigated the connection between phase response curves 
and mechanisms of brain function and pathology such as autoassociative memory jlOj and 
epileptic seizures |15j . The phase response curves of neurons in the neocortex have been 
experimentally determined [131 . 

1.2 First order phase reduction 

Before explaining our main results, which are collected in Theorem ll.il it is useful to briefly 
review the standard method of phase reduction. For a more complete discussion the reader 
is referred to [3, Chapter 10, and to [l]. 

It is easy to see why the derivatives of 8 are needed when studying weak perturbations 
of a stable oscillator described by i = F{x). Let such a perturbation be given by the new 
equation x = F{x) + eG{x, t), where e is a small positive number. Then 

^e(x(t]) = Ve(x(t)) ■i = u; + eVQ(x(t)) ■ G(x, t) (1) 
dt 

where V8 is the gradient of the phase function. Writing q']^ = V8(2:o(i)), where xo(i) is 
the point on the limit cycle on the same isochron as x{t)^ then 

^ = + eQ^P ■ G{xo{t),t) + eO(|x(0 - xo(i)l). (2) 
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One then proceeds by discarding the term eO{\x{t) — xo(t)|) and analyzing the resulting 
system. Thus implementing the phase reduction method in order 1 in \x{t) ^ XQ{t)\ requires 
finding the gradient of the phase function along the limit cycle of the unperturbed oscillator, 
denoted by Q^^'' in the last equation. 

One practical method for obtaining the gradient of O is by solving the equation: 

q[''> +DF\xo{t))Q[''> =0 (3) 

where the dot indicates time derivative and DF'^ {xo{t)) is the transpose of the Jacobian 
matrix of F evaluated along the limit cycle. This procedure was suggested by Malkin 
[71[TTJ[T2] and later by others independently [lH[3]. The reader is referred to [7], Chapter 9, 
for more details of the Malkin's theorem and [8], Chapter 10 for a historical note on phase 
reduction. One can find Q^^^ by numerically integrating the equation backwards in time for 
any initial condition satisfying ■ F{xo{0)) = 1, over an interval of time long enough to 
allow the solution to stabilize to a periodic one [17] . 

One limitation of the method as presented above is that it is only valid to first order. To 
develop higher order phase reduction, it is necessary to obtain higher order partial derivatives 
of O. However, to the best of our knowledge, an approach to higher order approximations 
of Q similar in spirit to the above due to Malkin for finding the gradient of Q has not been 
described so far in the literature. The goal of the present paper is to develop a numerical 
method to obtain partial derivatives of O to arbitrary order. 

Even within the framework of first order phase reduction (i.e. Equation ([2])), in certain 
situations, as when dealing with weak stochastic perturbations of oscillators, one may need 
to know the derivatives of along the limit cycle to orders greater than 1. For example, in 
[19| , a stochastic version of phase reduction is given using the 2nd order partial derivatives 
of O to obtain the mean and variance of the period of a limit cycle oscillator. They apply 
their results to the Stuart-Landau oscillator, for which an explicit form of can be obtained 
analytically. In general, however, an analytical form of is not available. Therefore, the 
method we present here is of particular interest in order two for studying the dynamics of a 
limit cycle oscillator subjected to stochastic perturbations. This will be discussed in detail 
in a future paper. 

The main result of this paper, which is described in Theorem 11.11 below . amoimts to a 
recursive procedure for finding the higher order derivatives of in which the first step is 
Malkin's method just described. 

1.3 The main result 

We need a few definitions first. Let v — (vi,...,Vn) be any vector in M" and / any 
differentiable real- valued function in xi, . . . , Xn- Given a multi- index J = (ji, . . . , j„), i.e., 
an n-dimensional vector with non-negative integer entries, we write = Vj^ . . . Vj^ and 
fj ~ D"^ f = Dl^ ...Dii'-f, where Di denotes partial derivative with respect to xt and 
~ Di . . . Di [k times). Let /^'^^ represent the symmetric /c-multilinear map on R" 
characterized (via polarization of polynomial maps) by 

\J\=k 

for all V G M". The sum is over all multi-indices J of order | J| = ji + • • • -I- jn = k. Here, 
and often later, we omit reference to the point x where the derivatives are taken. When 
necessary, this point is indicated as a sub-index; thus fx (f , w) is the bilinear map evaluated 
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at the vectors v,w regarded as tangent vectors at a; G M", where x is the point where the 
partial derivatives of / are calculated. 
We now define 

Qf^eJU. (4) 

for some fixed x e M", where 0t is the flow of F. Similarly, we define F^''''> , which is now a 
vector valued, symmetric, fc-multilinear map. (A convenient alternative description of G'-'^-' 
and F^'^^ , and more generally of the higher order derivative forms associated to tensor fields, 
will be given later in the paper.) In particular, F^^'^ is the linear map which to v G K" gives 
the directional derivative of F along v, i.e., F^^\v) = D.^jF = Tlij^i^i^ ■• where DyF is 
defined by this identity. 

Another general concept needed below is the symmetric composition of multilinear maps, 
which we define as follows. Let Q by a symmetric s-multilinear map on M" taking values in 
M, and H a symmetric fc-multilinear map on M" taking values in M". Then the symmetric 
composition of Q and H is the symmetric s + fc — f -multilinear map on K" , denoted Q & H 
and given by 

QQH{vi,...,Vs+k-l) = (gj^l_ 1)1 X! Q{H{v„(^i), . . . ,V„(^k)),Vai^k+l), ■ ■ ■ ,Va(s+k-l)) 

where the sum is over all permutations of the set {l,2,...,s + fc — I}. Finally, given a 
CO- vector Q on M" (a linear map from R" to M), we define the fc-multilinear map 

(Q ® • • • ® (5)(wi, . . . , Wfc) = Q(ui) . . . Q(ufc) (the fc-fold tensor product.) 

Theorem 1.1 Let x £ C, x{t) = (j)t{x), k > 1, and Qj'^' the k-multilinear map defined in 
(|4|) . Then, the following hold. 

1. Qt'^' satisfies the differential equation in ylj'''^ given by 

+ kAf^ Fif) = - E + 1) Q'r' F^:s' ' (5) 

where the right-hand side involves the Q\ for j < fc, and equals i/fc = 1. 

2. Let A he any k-multilinear map, N a positive integer, and At.N the solution to Equation 
(© for < t < NT such that A^t.n = A. Then there exists a T-periodic solution 
At such that At^N converges exponentially to At for 0<t<TasN—^oo. More 
precisely, there are constants C > and < A < 1 so that 

sup lAt^N-Atl^CX"". 

0<t<T 

3. Lf A^^'^ is any T -periodic solution of Equation then q'^^^ ^ A''^'^ -\- jiQ^^^ ®- ■ ■'SiQ')^ 
where 

= T''[Q^\f,, . . . , F,) - . . . , F,)]. 

4. The term QI^\Fx, . . . , Fx) has an a priori expansion as a linear combination of com- 
positions of the lower order terms Qq' and Fx^ for I < fc — 1. This expansion is 
described in section [KR 
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The equation for the standard (Malkin's) method of phase reductfon for obtaimng the 
gradient of Q is the equation in part 1 of the theorem when fc = 1: 

QP^+QP^°^S=0. (6) 

Furthermore, Q^^'' is the unique periodic solution such that q'^\Fx) = l/T. This unique 
solution can be found numerically, according to part 2, by the following procedure: Let a 
covector A (i.e., a linear map from R" to M) at p be a choice of initial condition for Equation[S] 
which is arbitrary except for the normalization A{Fp) = l/T. One then integrates Equation 
([5]) for t < (backward integration) until the solution stabilizes to a periodic (co)-vector- 
valued function on the limit cycle. Stabilization is assured to happen for sufhciently large 
\t\. This periodic function is the solution we want in order one. The general order case is 
then given recursively by the successive applications of the theorem. 

Before presenting the proof, we illustrate the use of the theorem with the Stuart-Landau 
oscillator. 



1.4 The Stuart-Landau oscillator: an illustration 

To illustrate the method, we focus attention on the case k = 2. From the general definition 
of symmetric composition introduced above we have that Qf^ F^"^^ and q'^P F^"^^ are 
given by 

QP' 0i^(i)(z;i,i;2) - \ (g(')(TW(i,i),z;2) + Qp)(F(i)(z;2),i;i)) 

We suppose that Q^^^ has already been obtained (say, by the standard method) and 
wish to find QJ . According to the main theorem, this second order term satisfies the 
non-homogeneous differential equation 

gf)-F2Qp)0f^W=-Q«0F(2), (7) 

where the F^^^ are evaluated at x{t) on the limit cycle. The equation can be solved as 
follows: Let be a solution to Equation [7] obtained by backward integration for an 
arbitrary initial condition. For large enough |t| this solution stabilizes to a periodic (tensor- 
valued) function along the limit cycle, which we still denote by ' . Then, by item 3 of 
Theorem II. 1[ 

Qf) = ylp' (q^;}\f,,F^)-A^;^\f,,F,)) Q« 0Q|^) (8) 

where ® Q^^H^i; ^2) = Q^P {vi)Q^p {v2)- We show later that the expansion referred to 
in item 4 of the theorem amounts in this case to Q^\Fx, F^) = —Q'l^\Fi^\Fx)) ■ Therefore, 

^(2) ^ ^(2) _ (qW(^(i)(^^)) + Af[F^^F^)) g« ® (9) 

is the solution we seek. 

We now recall the Stuart-Landau oscillator. (See [12] •) Define 
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We regard points of as column vectors: x — {xi^X2)\ where f indicates transpose. 
Let a, b be real constants and p{r) a smooth function of r > such that p{l) = 1 and 
= X > 0- Now define a vector field on by 



F{x) ^ AaX~ p{\x\)AbX. 



(10) 



Then it is easy to check that the differential equation x = F{x) has a hyperbolic stable limit 
cycle given by 5^ = {a; S 



\x\ ~ 1}. In fact, r = |a;| satisfies: 
r^r{l- p{r)) = -x(r - 1) + o(r - 1), 



showing that the limit cycle is approached for t > Q with Lyapunov exponent ~x- Special- 



izing to p{r) = r , then r — r{l 



is easily solved: 



r{t) = 1 + 



1 — r. 
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(11) 



where rg — ?'(0). With the coordinate change xi ~ cos{ip + blur) and X2 — sin((y5 + tilnr) 
we can write the solution to i; = F(x) explicitly in the new variables r, ip by setting 



if[t) = (^0 + (a - 



(12) 



as can be easily checked. Therefore, 9(a;) = Lpjln modulo integer translations. For this 
example, we can calculate the derivatives of Q(x) explicitly, and then compare them with 
the numerical values derived from Thcorcm ll.il 

Implicit differentiation gives the first and second order derivatives of along the limit 
cycle. We write 6^ = -DiO, Oy = DiDjQ, where Di is partial derivative in Xi. Then 



= (61,62) 



bxi + X2 xi — bx2 
27r|a;|2 ' 2tt\x\'^ 

Identifying 6^^' with the Hessian of 6, we can write 



(13) 



e(2) = 


611 


612 






621 


622 





2Tr\x\* 
xi 



27r|x|* 
2xiX2—b(x2—x'l) 



(14) 



The tensors Q^^"^ and Qf are similarly written. Let (^{t) ^ ipo + {a — b)t + 61nr(t). Then 



1(2) 



and 



(2) 



Qii(t) Qi2{t) 

Q2l{t) Q22{i) 



bcos(2C(t))+sin(2C(t)) 

27rr2(t) 
- cos(2C(t))+fcsin(2C(t)) 

2ii-r^(t) 



- cos(2C(t))+bsin(2C(t)) 
27rr2(t) 
fccos(2C(t))+sin(2C(t)) 
2^7^7*1 



For any vectors v,vi,V2 £ R^, 

F = AaX - \x\'^Ai,x 
F^'^\v) = Aav - 2x ■ vAi,x - \x\^Ai,v 

F''^^Vi,V2) — —2x ■ ViAbV2 - 2x ■ U2^bWl - 2ui • V2AbX. 



(15) 



(16) 



(17) 
(18) 
(19) 
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Let the components of these tensors relative to the standard basis 61,62 of be denoted 
as follows: 



Flix) 
Ff{x) 



FUx) 



where the entries are obtained from Equations [iTl [THl and [191 For example, from Equation 
[TOlit follows that F^^{x) = F^^ix) = 62 • i^^^H^i, 62) = -2(a;i + hx2). The other entries are: 



Fl{x) 


F^{x) 




1 — 2xi{xi — 


bx2) - |X|2 - 


-6- 


2X2(X1 - 


- 6x2) + 6|xp 


FHx) 


Fi{x)_ 




b — 2xi{hxi H 


-a;2)-6|x|2 


1 - 


2x2(6x1 


'rX2) - |xp 


Fhix) 


Fhix) 




— 6x1 + 26x2 


26xi — 2x2 








F^iix) 


FUx)_ 




26x1 — 2x2 


—2x1 + 66x2 








Fhix) 


FUx) 




—66x1 — 2X2 


— 2xi — 26x2 








Flix) 


FUx) 




— 2xi — 26x2 


— 26xi — 6x2 
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Time (sec) 

Figure 1: Determination of Q'^' for the Stuart-Landau oscillator with a = 2, 6 = L Each compo- 
nent of the tensors is plotted as a function of time. Numerical integration is done backward in time 

using the Euler method with dt = \Q . Solid line: QJ determined from At,NT', dashed line: Ql 

(2) 

found analytically. To clearly show convergence of the numerically obtained Ql to the analytical 
solution. At, NT was not plotted and the period of the simulation was set at only 0.8T. 



We can now write the differential equations for Qi and Q.y (Equations [6] and [71 all 
summations are over s ~ 1,2): 
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Q.^Yl Q^P^ (20) 

S 



Q^^ = T.^f:Qs, + f;Qs^) - QsF^ . (21) 

s s 

The first equation is tlie one used in tlic standard Malkin's approach. Once it is solved 
by the aheady indicated procedure, its solution enters as the non-homogeneous term for the 
second equation. 

The result of the numerical calculation is shown in Figure [1] The initial condition for 
^t.NT was chosen from a set of random numbers and numerical integration of Equations 
(|2ip was done backward in time (data not shown) . Then, Qi was obtained using Equation 

([8]) (solid line in Figure [1]). Convergence of Qp', obtained numerically, to the analytical 
solution (dashed line in Figured]) is clearly observed. 

2 Proof of the main theorem 

The subsequent sections are dedicated to proving Theorem ll.il It is convenient, for reasons 
we hope will become apparent along the way, to use a more geometric and coordinate- 
free language for manipulating tensors, even though all calculations are done in K". This 
preparatory material on tensor calculus will take a few pages to develop, so it may be 
appropriate to provide an overview of the proof in simpler terms first. 

The basic idea for proving part 1 of Theorem 11.11 is the following. An immediate con- 
sequence of the definition of the phase function is that G(0t(7(s))) ^ uit + 6(7(3)) for any 
differentiable curve 7(5). Therefore, .A^e(^j(^(s))) = for fc > 1. This gives an ODE 
that the kth derivative of the phase function must satisfy. Thus the proof of part 1 simply 
amounts to successive applications of elementary differentiation rules, although handling the 
multitude of terms that arise quickly becomes a challenge. In fact, this is precisely the kind 
of situation that often calls for the use of Faa di Bruno type formulas. We take here a differ- 
ent approach, which obtains the differential equation recursively and uses a coordinate-free 
language to facilitate the manipulation of terms. 

The proof also involves a somewhat surprising cancellation of terms that renders the 
result more simple than might be expected at first. To better understand this point, we 
sketch the derivation of the differential equation in dimension one. Although the calculation 
is considerably more complicated in the general case, the key formal aspects of the derivation 
are already present for n — 1. We use subscripts to indicate the number of derivatives, so 
fk{x) is the fcth derivative of f{x), whereas f'{x) is also used for the first derivative. Let 
7(5) = xo +VS, where xq and v are fixed. Let Ut{s) = ^4(7(5)), so u[{0) = (/)J(.to)w. Also set 
qI''\s) = Ofc(ut(s)). Recall that we wish to find a differential equation for q[''^ = q[''\o). 
The differential equation associated to F is ut{s) ~ F{ut{s)). Then for k = 1, 

= Jt^,^iMim = {Qi'\'-^) + Q['\s)F,iu,{s))} 4ilis))v 

for arbitrary v. An exchange of the order of derivatives in < an s was used in the second term 
of the right-hand side of the equation. Thus we conclude that Q[^\s)+Fi{ut{s))Q[^\s) — 0, 
from which the claim for fc = 1 follows by taking s = 0. 
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We now suppose that 



for 1 < J < fc and wish to show that this also holds for j — k + 1. Taking one more derivative 
in s and exchanging derivatives in t and s gives 

" = i {l®'^-^''*^'^^ + ^ QFKu*(,s))e,+i_,K(s))| 

= {Q[''+'\s)+Q[''+'\s)F,iu,{s))}u[is) 
1=1 ^ ^ 

Thus the desired equation will hold for j = fc + 1 if we can show that 

1=1 ^ 

This indeed holds and is easily checked for n = 1 by using the well-known relation 

The n > 2 case is considerably more involved but follows a similar line of proof. 

Our use of tensor calculus becomes more essential in the proof of parts 2 and 3 of the 
main theorem, which rely on a study of the stability properties of the tensor equation given 
in part 1. Part 4 boils down to an enumeration of certain combinations of tensors that can 
be represented by formal linear combinations of rooted trees. This is discussed later under 
the heading "forest expansion." 



2.1 Higher order differentials of tensor fields 

In this and the next two sections wc develop some background material on differential 
calculus of tensor fields that will make our calculations of Taylor expansions of functions 
and vector fields more tractable. We suppose that the standard concepts (see, for example, 
[S]), such as tangent and cotangent bundles. Lie derivatives, covariant differentiation, etc, 
are known, but recall some of the definitions for the purpose of setting up notation. 

Let y be a finite dimensional vector space and V* its dual space. The space of tensors 
of type (r, s) is 

= v ®---(^V (^V* ®---®V\ 

in which there are r copies of V and s copies of V* . We think of an element of 
as an s-multilinear function taking values in the set of (r, 0)-tensors. We find it useful to 
represent a tensor T of type (r, s) as a diagram consisting of a box with s (contravariant) 
lower legs and r (covariant) upper legs, and think of the lower legs as places where the vector 
arguments are plugged in. See for example, Figure [H 
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A tensor field of type (r, s) on M" is a function that associates to each x in M" an element 
of Vx^'''\ where = T^-M" is the tangent space at x. (This tangent space is, of course, 
naturaUy identified with M", but observing the distinction will help keep track of where 
derivatives are evaluated.) 

All the tensor fields (functions, vector fields, etc.) that are considered below are smooth. 
The result of evaluating an (r, s) tensor T on s vectors Xi, . . . ,Xs is the (r, 0) tensor denoted 

r(ui, . . . , Vs). 

Let D be the standard covariant derivative in R". If X is a vector field and u is a tangent 
vector at some point, then DyX is the standard derivative of X along v. The derivative of an 
ordinary function / along v will be written as vf = Dyf. We can then extend the definition 
Dy to general tensors by imposing the condition that the product rule for differentiation 
holds for the tensor product CE) as well as for the pairing operation of V and V* . This implies 
that ii Xi, . . . , Xs are vector fields and v G T^M" , then 

s 

{D,,T){Xi, ...,Xs)^ D,T{Xi, ...,X,) - ^T(Xi, . . .,DyX,, . . . ,X,). (22) 

When T has type (0, s), T{Xi, . . . , Xg) is an ordinary function, and we write vT{Xi, . . . , Xg) 
instead of DyT{Xi, . . . , Xs). The result of the differentiation of T is an (r, s + 1) tensor field 
DT. 

Similarly, D^T is the (r, s + fc) tensor field obtained from T by applying D fc-times. 
We call D'^T the fc-th order differential of T. Regarding TC^') = D^T as a fc-multilinear 
map taking values in the space of (r, s) tensors, then it can also be defined inductively: 
r(i)(u) ^ DyT and T('=)(«i, . . . ,Vk) = {Dy,T^''-^^){v2, . . . , ffc) for fc > 2. 

We are dealing with the standard covariant differentiation in M", so D is symmetric (i.e., 
its torsion tensor is 0) and ffat (i.e., its curvature tensor is 0): given vector fields Xi, X2, X3, 
then DX1X2 - DX2X1 = [Xi,X2] and Dx.Dx^Xs - DX2DX1X3 = D[Xi^X2]X3, where 
[Xi, X2] is the Lie bracket. We also say that an (r, s)-tensor field S is symmetric if for each 

permutation a of the set {1, . . . and each x, S'i'^'' (wcr(i) 7 • ■ • iVa{k)) = si''\vi, . . .,Vk) for 
all vectors wi, . . . , w^, £ Tj-R". For simplicity, we are not indicating the s vector arguments 
of S. 

Proposition 2.1 Let T he an {r, s)-tensor field in R" and T^'^'> the k-th order differential 
of T . Then T^''^ is symmetric. 

Proof. We indicate the proof for k = 2. The general case follows by induction using the 
same argument. Let Xi,X2 be vector fields that agree with vi,V2 at x. If follows from the 
definitions that 

T^^\vi,V2) ^ Dx.Dx^T - Dd^^x.T 
where the right-hand side is evaluated at x. Therefore, 

Ti^) (vi ,V2)- TP {v2,v,)^ [Dx, , Dx2]T - ^ = 0, 

where the last equality is due to the fact that D is ffat. □ 

As an example, let T be an ordinary function, denoted /. Let Xi, . . . , X„ be the coor- 
dinate vector fields in R" : Xj ~ g|- . Then it is easy to check that 

Z^'"' i^ji 1- ■ ■ y ) = ^ji ■ ■ ■ /i 

which is the A:-th order partial derivative of / with respect to Xj-^ , ■ . ■ , Xj^. . 
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2.2 Symmetric composition of tensors 

Let Sk represent the group of permutations of the set {1, . . . , k}. Given a tensor T of type 
{r,k), we define its symmetrization T as the (r, fc)-tensor derived from T by symmetrizing 
its k contravariant legs: 

T{vi,...,Vk) = X! • ■ • ^■''fffc)- 

If we need to symmetrize only a subset of the contravariant legs, this will be indicated in 
some explicit fashion. For example, if T is an (r, k + /)-tensor, we separate by a semicolon 
the vector arguments that will not be symmetrized and place them last in order of insertion: 
T{vi, . . . , Vk]Wi, . . . , wi). Diagrammatically, the input legs taking the arguments wi, . . . ,wi 
could be shown, for example, on either side of the tensor box, although we will simply omit 
them in our diagrams. A thick line crossing the other contravariant legs is added to indicate 
symmetrization. 



T 



Figure 2: Diagram for the symmetrization of a (r, fc)-tensor. The symmetrization of the contravari- 
ant legs is indicated by a thick bar. 

It is convenient to introduce a binary operation on tensors, which we call symmetric 
composition. Let Q be an (r, s)-tensor Q and H a fc)-tensor, where I < s. We define the 
symmetric {r,s + k — Z)-tcnsor Q Q H as follows: 

QQH{vi,...,Vs+k-l) = tg^l_ ni QiH{Vail), ■ ■ ■ ,Va{k)),Va{k+l), ■ ■ ■ 

We have assumed here for simplicity that all the contravariant legs are involved in the 
symmetrization. The more general case is defined similarly. 

It is clear from the diagram that Q Q H would still make sense if ? > s, in which case 
the result of the operation is a symmetric (r + / — s, A:)-tensor. Below we mainly need the 
case where / = 1. We also note that the input legs to which symmetrization is typically 
applied later in the paper are created from multiple applications of covariant differentiation 
to a tensor. 

Proposition 2.2 The symmetric composition of Q and H satisfies the product rule: 

D^iQ QH)^ {D,Q) QH + QQ D,H. 

Proof. The proof is a tedious but straightforward application of Equation [22] and the easily 
proved fact that the {k, fc)-tensor P : vi (E) ■ ■ ■ (g) Vk — > "(^(7(1) <8) • • • ® ^^(fe) satisfies DP = 0, 
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Q 



H 



Figure 3: Symmetric composition Q Q H oi Q and H . 



where a is any permutation. 



□ 



When Q and H are themselves symmetric, Q Q H has the following useful description. 
Let C(fc, I) represent the collection of all /-subsets of {1, ... , k} (i.e., subsets with / elements). 
The cardinality of C(fc, I) is the binomial coefficient (^) = k\/{k — l)Ul. Let ir : Sk —* C(fc, I) 
denote the map that associates to each permutation a the set J ~ {o'(l), . . . , cr{l)}, and write 
Sj = TT~^{J). Then each Sj, J G C{k, I), has cardinality U{k — 1)1 and Sk is partitioned into 
the disjoint union of the Sj. The following notation will be used: let wi, . . . , Wfe be a subset 
of TJeR". H a symmetric (m, /)-tensor at x, and J an /-subset of {1, ... , k}. Then it makes 



for J ~ {ji, . . . ,ji}. If Q is a symmetric (r, s)-tensor. 



sense to write H{J) 

where s > m and s ~ m = k — I, we similarly write Q{H{J), J"^). where J"^ indicates the 
complement of J. (More generally, Q and H may contain additional contravariant legs, 
with respect to which the tensors are not necessarily symmetric, that arc not included in 
the operation of symmetric composition.) 




Q 



H 



Q 



H 



B 



C 



Figure 4: Tensors A, B, C satisfying the identity {i^{)A + ( 7 )B = of Proposition lOl 



If i; G TxMJ^ and H is an (/, fc)-tensor, let iyH denote the {l,k ~ l)-tensor such that 
{iyH){v2, ...,Vk) = H{v,V2, ■ ■ . ,Vk). 
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Proposition 2.3 Let H be an (u, l+l)-tensor, Q an (m, k—l — l+u)-tensor, both symmetric 
in their lower legs, and v G T^K'^ . Then 

Proof. Let vi — V, V2, - ■■ ,Vk be elements ol TxWK The key point to observe is that the 
sum ol Q{H{J), J'^) over all J £ G{k, ^ + 1) is equal to the sum over all J S C(fc, I + 1) that 
contain 1 plus the sum over all J that do not. □ 



2.3 Identities for the Lie derivative 

The infinitesimal action ol the flow 0t on tensors is given by the Lie derivative with respect 
to the vector field F. We register here lor later use some uselul formulas involving the Lie 
and covariant derivatives. We recall that the Lie derivative ol tensor fields is defined just as 
we did for the covariant derivative, except that, on a vector field X, it is given by the Lie 
bracket LpX = [F,X], and on functions Lpf = Ff. Another way to describe UpT is as 
the time derivative ol the tensor defined by "convecting" T along the flow of F. 

Let now be a smooth vector field (thus F is a (1, 0)-tensor) and set F^''^ ~ D^F. Thus 
F'*^' is a (1, fc)-tensor. For the moment, we only make use ol F^^\ 

The covariant and Lie derivatives of tensors with respect to a vector field F are related 
via an "algebraic derivative" with respect to F^^' . The latter operation is defined as follows. 
Let i? be a field ol linear maps and define Ab by the properties: (1) AbJ = on lunctions; 
(2) Abv = Bv on vectors; (3) AbO- = —aoB on covectors; (4) Ab satisfies the product rule 
with respect to tensor product. It can be checked that these properties uniquely determine 
a linear map on tensors. In particular, if Q is a (0, A;)-tensor, 

fc 

(■AbQ)(ui, . . . ,wfc) = - ^Q{vi, . . . , Bvj, . . . ,vk). 
i=i 

From the properties stated in Proposition 12.41 below we obtain that il Q is a (0, A:)-tensor 
field, then for any point x and tangent vectors ui, . . . , Ufc at x, 

fc 

{^FQ)xivi, ...,Vk) = {DpQ)x{vi, . . . , Wfc) + ^ Qx{vi, ■ . .,F^'^\vj), . . . ,Wfc), 

and il moreover Q is symmetric, then HpQ — DpQ + kQ Q F^^\ The next proposition sum- 
marizes some relations among the various derivative operations. Here [■, ■] is the standard 
commutator ol operators. Also note that, il we regard F as generating the time evolution 
ol a system, then with a slight abuse ol notation, DpQ = Q, where the dot over Q means 
time derivative. Below, F^'+^)(Xi, . . . ,Xi,-) is the field ol linear maps defined by inserting 
a vector into the last slot. 

Proposition 2.4 Let F, X , Xi, . . . ,Xi be vector fields and Q an (r, k)-tensor field. Then 

1. DpQ- LpQ^ApwQ 

2. [Dx-,Ap(i)]Q=Ap(2)(^x,-)Q 

3. [Dxi,Ap(i)(^x2....,Xi,-)] = ■^F^'+^^Xu...,Xl,■)+J2'-^^ F('){X2,...,DxiXj,...,Xi,-) 
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4- i^F, Dx]Q — D[p x]Q + ■A.p(2)(^x.-)Q 

5. {ZpQ^^^){X; .) = DxiipQ + ^fp)(x,)Q- 

Proof. It is easily checked that Dp — Hp has the properties defining Ap(i), so 1 is a con- 
sequence of uniqueness of Ap{i). A similar verification also proves 2. Property 3 is shown 
by induction and the same argument used for 1 and 2 based on uniqueness of the algebraic 
derivative. Properties 4 and 5 can be proved using 1 and 2 by a tedious but straightforward 
algebraic manipulation. When deriving these properties, it should be born in mind that D 
is torsion-free and flat. □ 

Property 5 of the above proposition gives a way of finding LpQ^^^ recursively if HpQ is 
known. We illustrate this with a formula for HpQ when Q = f'^^^ and / is a function. This 
is the case we need to consider in extending Malkin's method. The key point to notice is 
that the formula expresses Lpf'^^'> in terms of the lower order tensors /'■'■', / = 1, . . . , /s — 1. 

Proposition 2.5 Let F he a smooth vector field and f a smooth function. Suppose that 
F f = g, for some smooth function g. Then Lpf^^'' = g^^^ and, for k >2, 



1=2 ^ ^ 



Proof. The proof is by induction. Using Property 5 of ProDOsition l2.4l one immediately gets 

^Ff^^'> = g^^^ and 

[Zpf^'^) («1,«2) = {D,,Zpf^^^) («2)+ («2) 

= {g^''^-f^''^&F(^^)iv,,v,). 

So we suppose that the equation holds for /c > 2 and wish to obtain it for fc + 1 . First note 
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that 



= g('=+i)(«i,...,T;fc+i) 



1=2 

k 



;=2 ^ ^ 



fe-i 



fc-i 



1=1 



+ fc(/('=)0i?,,FW) (?;2,...,t'fe+i). 
Using the tensor identity of Proposition 12.31 the above simpUfies to 

1=2 ^ ^ 



+ fc(/W0i?,,F(l)) (t.2,...,«fe+l). 
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By Property 5 of Proposition [^131 

k+1 

;=2 ^ ^ 



fe+i 



i=2 

fe+i 



( — 9 \ / 



1=2 

which is the claimed formula for fc + 1. □ 

The case of main interest is / = O, for which F6 ~ l/T is constant. Thus g^''^ = for 
fc > 1 and we obtain 

LfQ' = -^2(1) ©*'+'"'^ f'^'^- (23) 



1=2 



2.4 Proof of Theorem [TT] 

We turn now to the proof of Theorem 11.11 Part (1) of the theorem, which gives the dif- 
ferential equation satisfied by O'^'^-' is restated in the next proposition. Recall the notation: 
Qj*^' = ©^'^l^'), where (j)t is the flow of F and x is arbitrary. Indicating the time derivative 

of QC^) by g("\ we have Q^'^ = (i^j. e(^-))^,(,). 

Proposition 2.6 We assume the notation of Section \1.3[ Let Q'^^^ he as just defined. Then 



(=1 



Proof. According to Proposition!^ DpQ^^^ - Api^Q^^'^ = LpO^'^l On the other hand, 

since the directional derivative FQ is constant, ZLj^B^'^-' is 
given by Equation l23l We get the claimed formula by finally rewriting the resulting equation 
in terms of Qj'"' . □ 

In what follows, let Wt be the tangent space to the isochron at (j>t{x), x £ C. Vectors 
in Wt, by assumption, contract exponentially under the flow; i.e., \d(j)sv\ < CA*|v| for all 
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V G Wt and positive constants C and A < 1, where s > 0. If v is parallel to F{(j)t{x)) 
then \d(l)sv\ is bounded above as well as away from 0. We similarly need to know the decay 
properties of tensors of type (0, k) under the flow. The natural push-forward action of 0s 
on a tensor- valued function, Tf , of type (0, k) along C is defined by 

where ui, . . . ,Uk are vectors at y. This applies, in particular, to a tensor field r defined in 
a neighborhood of C, in which case r-t = T^t{x) is a T-periodic tensor- valued function of t. 
The function is fiow- invariant if (jjg ■ Tt = Tt+s ■ 

Let Mr represent the one-dimensional space over M spanned by a tensor t. Since the 
family of isochrons and the vector field F are invariant under the flow, the tangent space to 
M" decomposes invariantly as a direct sum 

where Ft = F{4)t{x)). Let Wj* be the subspace of T^^^^^M" consisting of covectors that 
vanish on Ft. It is not difficult to check that 

is also a flow-invariant decomposition. By general tensor algebra, one also obtains a flow- 
invariant decomposition of the space of (0, fc)-tensors as a direct sum of subspaccs of the 
form 

Vk,M,t = KQ^^' ®---® rqI^^ ®w; ®---® Wt*, 

in which there are ki copies of RQI^'' and k2 copies of Wt , where ki + k2 = k and < fci < fc. 

The Euclidean norm on vectors, extends in natural ways to norms on tensors of 
any kind. We use the same notation, |r|, for the norm of a tensor r of general type. We 
refer the reader to texts on multilinear algebra or differential geometry for how this can be 
defined, although it is not necessary for what we do below to have any explicit description 
in mind, and the form of the theorem does not depend on a particular choice of norm. The 
main property we use below is that if A is a (0, fc)-tensor and ui, . . . are vectors then 

Lemma 2.1 Let A he a tensor in Vk-^^k2,t- Then there exists a constant Ca such that 
10s -A] < Ca A'=2l"l /or a?; s < 0. 

Proof. A tensor in Vk^,k2,t has the form q[^^ ® ■ ■ ■ ® Q^^^ ® A where A belongs to the 
fc2-fold tensor power of Wf. Thus |0s • A| < K\(f)s ■ A\, where K is the fcith power of the 
supremum of IQI^^I over the period < s < T. Since A vanishes whenever any one of its 
arguments is parallel to Ft , the norm of A can be bounded above by the supremum of 
\A{{d(l)s)~^ui, . . . , {d4>s)^^ui)\ over all Uj 6 Wt+s of norm at most 1. Now, 

which proves the assertion of the lemma. □ 
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Proposition 2.7 Let a; G C and A[^^ a (0, k)-tensor at (jit^x) depending differentiably on t. 
Suppose that A^^^ is a solution of the homogeneous equation 



A 



Then for some positive constant K and for \ < \ as above, there exists c such that 



A 



(fc) 



(1) 



< 



for all t < 0, where the tensor product on the left-hand side contains k terms. In particular, 



if A^''^ is T -periodic, then A^^^^ 



Proof. Since Fx — along the limit cycle (by the usual identification of a vector field 
as a derivation), At equals Dp A, where the latter has to be interpreted in terms of an 
arbitrary extension of A on a neighborhood of the limit cycle. Now, according to part (1) 
of Proposition 



LpA 



(fc) 



a: 



■ Ap{i) A 



(fc) 



i(fc) + fc^(fc) pii) ^ 



on C. Therefore, the proposition amounts to the assertion that (1) the space of periodic, 
flow-invariant tensor-valued functions of type (0, fc) on C is one-dimensional, spanned by the 
fcth tensor power of Q[^\ and (2) if A'f ^ is flow-invariant but not necessarily periodic, the 
stated inequality holds, i.e., the components of a solution transverse to the one-dimensional 
space of periodic solutions must contract exponentially to 0. 



Thus suppose first that A], is flow-invariant and T-pcriodic. Let A 



(fcl,fc-2) 



be the com- 



ponent of j4j'^^ in Vfcj,fc2,f. In particular, A'f^'^'' 



ctQ\ 



(1) 



® Q^^' . Since the subspaces 



(fel,fe2) 



^fci,fc2,t flow-invariant and periodic, each component A. 
periodic. In particular, q is constant (independent of t) and A\ 



is also flow-invariant and 



(fcl,fc2) 



by Lemma [2.11 This shows claim (1). To verify claim (2) we suppose that A\ is flow 



whenever fc2 7^ 0, 

(fc) 



invariant but not necessarily periodic. Thus A^^^s ^ 
decomposition 



is flow- invariant, so ct+s = ct for all s and 



1(1) 



(/)s • A^^^ . Now each component in the 



< 

fc2#0 
fc2#0 



(fel,fe2) 



■A 



(kiM) 



for some constant K, where we have used again Lemma [2.11 This concludes the proof of 
the proposition. □ 



We can now prove parts (2) and (3) of Theorem ll.il For fc = 1, the equation reduces to 



Q«oFW=0. 
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This equation can be solved numerically by the standard (first order) approach, which 
essentially amounts to Proposition 12.71 in the special case k — 1: Fix x £ C and let A be a 
choice of initial value, arbitrary except for the condition A{Fx) — 1/T. One then integrates 
the first order equation for t < for large enough values oi \t\ until the solution stabilizes 
to a periodic matrix-valued function over the limit cycle. That function is the sought after 
solution for Q[^\ For larger values of k, we regard the differential equation in Proposition 
12.61 as non-homogeneous: 

where the right-hand side is assumed to have already been obtained in the previous steps. 
We proceed by induction. Suppose that we have found q[^\...,q['^ and wish to 

(k) 

obtain Q). ' . Set 

Bl'' = -j:(il,)Qt''QF^'^^^ (24) 
1=1 ^ ^ 

and integrate the equation 

(k) 

for t < 0, starting with an arbitrary initial condition, until stabilizes to a periodic 
function of t. Stabilization must occur at the exponential rate given by Proposition 12.71 
since the difference A^'^'' — Q^^^ is a solution of the homogeneous equation of Proposition 
12.71 This periodic solution is still denoted by A^^^ (in particular, A^^^ is the value at x after 
stabilization) . The true solution we seek, , differs from by a periodic homogeneous 
solution so there exists a constant c such that 

Since Q^^'\f) = l/T, the constant c has the form 

c = T^[qW (F„ . . . , f^,) - 4'=) (F„ . . . , F,)]. 

It remains to argue that Q^q\Fx, . . . ,Fx) can be expressed in terms of the lower order 
(5q\ F, and its derivatives at x. Recall that Q^^^ = ©^'^l^.-j- The following general fact holds 
for Q^y^ (for y not necessarily on the limit cycle). We set = F. 

Proposition 2.8 There exists for each k an algebraic function that gives Qy^\Fy, . . . , Fy) 

in terms of the &y\ for I = 1, . . . , k — 1, and the Fy\ for j = 0, . . . , fc — 1, and y in a 
neighborhood of C. 

The precise meaning of this proposition, and an algorithm for obtaining the indicated 
algebraic function, are explained and illustrated in the next section. 

2.5 Forest expansion of F'^Q 

We wish to expand 8'^'^-'(F, . . . ,F), at any given point, in terms of the O'-'' and F'^^^ for 
J < A:, and thus find the algebraic function indicated in Proposition l2.8l This can naturally 
be done by successive applications of the chain rule, applied to the expression F'^Q = 0, 
k>2. (Recall that FQ = O^^^F) = 1/T.) What is needed is a method to conveniently deal 
with the combinatorial task of enumerating the terms of this expansion. This can be done 
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by enumerating rooted trees with k edges. We describe here this forest expansion method 
and illustrate it with a few examples. 

Finding the general form of the algebraic function claimed in Proposition 12.51 for an 
arbitrary k is a complicated combinatorial problem, and amounts to a type of Faa di Bruno 
formula, which we do not attempt to describe here. We are content with giving the forest 
expansin algorithm and applying it to small values of k. 

Figure [5] explains how to represent nested evaluations of tensors by rooted trees. The 
number of edges of the rooted tree diagram is the total order of differentiation, so all 
diagrams associated to F'^Q will contain k edges. The degree of the root vertex is the 
order of differentiation of 8 (this is 3 in the example of Figure [5]) , and all the other vertices 
represent a where j + 1 is the vertex degree. Thus to each leaf (i.e., a vertex with 

no descendants) is attached an F, and to each non-root vertex is recursively attached a 
vector as follows: Starting from the leaves (associated to copies of F) one moves down 
one step to the parent vertices (associated to copies of F''^\ if j + 1 is the degree of a 
parent vertex) and compute F'-^^F, . . . ,F). This vector is now attached to each of those 
second-to-last generation vertices. These new vectors in turn are evaluated into the tensors 
attached to their parent vertices. We continue this process until the vectors attached to the 
first generation vertices (the ones connected to the root by an edge) are evaluated into 0'-'^ 
where I is the degree of the root vertex. 




Figure 5: A rooted tree diagram representing Q^''^\F^^\F^''^\F, F, F), F'^>(F)), F, F), which is one 
term in the expression of F^Q. The tree itself will be represented by T3(r2(r3, ti), tq, ro). 

The expression 9^'^-'(F, . . . ,F) itself corresponds to a tree having a root of order k and 
k leaves attached to it. We denote this tree by Tk, k > 0, where tq consists of a vertex with 
edges, i.e., a leaf. Other trees are obtained by nesting trees of type tj. Thus each r, can 
take j arguments, each of which is a tree of the same kind (for possibly different j). For 
example, t; (t^j , . . . , ) represents a tree that consists of a root vertex of degree / and at 
the non-root vertex of each of the I edges is appended the tree r,^ so that the root vertex 
of the latter is identified with the non-root vertex of the former. 

It is clear that F^Q can in general be represented by a forest of rooted trees with k 
edges, each tree being assigned some multiplicity. We will see shortly how the multiplicities 
are determined by counting the ways a tree is derived from other trees with fc — 1 edges. 
Figure [2T5l shows the forest diagram representation of F^Q, F'^Q and F'^Q. 

We give a few examples of the forest expansion before showing the general method. 
By the basic rules of covariant differentiation of tensors, we obtain FF& = i^0^^^(F) = 
9(2) (F,F) + e(i)(i^(i)(i^)). Therefore, 

e(2)(F,F) = -e(i)(i^(i)(^))- 
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For fc ~ 3, wc have: 

FFFQ = i^[e(2) (F, F) + e^i) (F))] 

so that 

-e''^\F,F,F) = 3e(2)(F^i)(F),F) + e(i)(F(2)(F,F)) + e(^)(F(i)(F(i)(F))). 

The expansion of — 0'*'(F, F, F, F) is shown diagrammaticaUy in Figure [6l 




Figure 6: The forest representation of -O^^^F, ■ ■ ■ , F) for fc = 2, 3, 4. 



It should now be apparent that the algebraic expression giving — 8*^'^^(F, . . . , F) as a 
function of the lower order terms O*^-'-' and the F^^\ as claimed in Proposition l2.5l is precisely 
the forest expansion of —&^^'>{F, . . . , F). It is also clear that all rooted trees with k edges 
(except Tfc), up to isomorphism, appear as a term in the forest expansion of Tfc, for a given k. 
What is needed then is a more formal description of how to determine the integer coefficients 
of the expansion. (Notice our slight abuse of language in referring to the forest expansion 
of F^e or of -e'^'^^F, . . . , F) as the same thing. Of course, the expansion of the former 
contains one extra term, which is (minus) the latter.) 

Let T(fc) denote the set of (isomorphism classes of) rooted trees with k edges, k > 0, 
and denote by m : 7{k) — > N U {0} the multiplicity function^ which assigns to each tree 
its coefficient in the forest expansion of F'^Q. As already defined, each tree gives rise to a 
number: to the root vertex we associate 0*-'^ where I is the vertex degree; to each of the 
other vertices we associate F^^\ where j is the degree of the respective vertex and F^^^ = F; 
the tensors are then evaluated as prescribed by the tree so that each vector attached to a 
vertex is an argument of the tensor attached to the parent vertex. The result of this nested 
evaluation of tensors for a given T € 'J'(fc) will be written T{Q, F). Therefore, 

F'=e= J2 MT)Tp{e,F). 

TeT(fe) 

Thus, the forest expansion of F'^0 requires an enumeration of all rooted trees with k edges, 
and the determination of the multiplicities m{T). 

A few more definitions are needed before identifying m. A tree F' £ T{k) is said to grow 
into T £ T(A: + 1) if F can be obtained from T' by adding (grafting) one terminal edge to 
any vertex of F'. Conversely, F is pruned down to F' if T' results by eliminating a terminal 
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edge from T. Let Cfe be the N-module consisting of linear combinations over N of elements 
of T(fc). The dual of T G T(fc) wiU be written as T*, so that T*{Y,s(^7(k) "^s^) = «t- 

The pruning map V : 'J(fc + 1) — > 6^ is defined as follows: For each T E 7{k + 1) we 
set J'(T) to be the sum of all distinct trees in T(fc) which can grow into T . The grafting 
map S : Cfc — > C^+i is defined on a tree T by summing all trees that can be pruned down 
to T, now counting repetitions (i.e., each tree is counted as many times as it appears in the 
process of grafting an edge at the different vertices), then multiplying the result by m(T). 
Now extend S to Cfe by linearity. It can be shown that the multiplicity function has the 
form: 

m{T)^T*{c,{y{T))). 

This is based simply on translating into diagrams the rule: If ^ = F'^^\Ai, ...,As), then 
Dp A = + Yf.^^F^''\Ai,...,DFAj,...,As). The first term on the 

right corresponds to adding an edge to the root of the tree associated to A; the other terms 
on the right correspond to moving up one step along one of the root edges of A and repeating 
the operation. 

We illustrate the procedure for finding multiplicities with the example of Tfe(ri, tq, . . . , tq) 
shown in Figure [71 We use the shorter notation Tk,i for this tree. 




Figure 7: The rooted tree Tk,\ ~ Tk{Ti, ro, . . . , tq). The pruning of the single edge at distance 2 
from the root gives Tk, and the pruning of any of the other edges gives Tk-i,i = rfc_i(ri, ro, . . . , to). 

We first calculate m{Tk). Clearly. J'(Tfe) = Tk-i- The grafting map gives: 

Sinn)) ^ s{Tk-i) 

= m(Tfc_i)((fc - l)rfe_i,i + Tk) 
= m{Tk-i)Tk + 

Therefore, m{Tk) = (S(J')(ta;)) = m{Tk-i)- It is clear that ?7i(ri) = 1, so m{Tk) = 1 for all 
k. 

We apply the same argument to rk,i. First, T'(Tfc.i) = + t^-i^i. Now, 

s(a'Ki)) = s(T,) + g(T,_i,i) 

= m{Tk){kTk,i + Tk+i) + m{Tk-i,i){Tus + ■ ■ ■) 
= {km{Tk) + m{Tk-i.i))Tk,i + ■■■ 
= {k + miTk-i,i))Tks + ■■■ 

Denoting by ruk+i the multiplicity of Tfe^i, we obtain ruk+i = ruk + fc, TO2 = 1. A simple 
induction gives the result: m{Tks) = k{k + l)/2. 

It would be useful to derive general properties of the multiplicity map that can help to 
evaluate the forest expansion of -e^*-') {F,...,F). For example, if T is any rooted tree, then 

m(r) = m(ri(---Ti(r) •••)). 
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